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Synchronization is optimal in non-diagonalizable networks 
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We consider the problem of maximizing the synchronizability of oscillator networks by assigning weights 
and directions to the links of a given interaction topology. We first extend the well-known master stability 
formalism to the case of non-diagonalizable networks. We then show that, unless some oscillator is connected 
to all the others, networks of maximum synchronizability are necessarily non-diagonalizable and can always be 
obtained by imposing unidirectional information flow with normalized input strengths. The extension makes 
the formalism applicable to all possible network structures, while the maximization results provide insights into 
hierarchical structures observed in complex networks in which synchronization plays a significant role. 
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Under extensive study in recent years is how the collective 
dynamics of a complex network is influenced by the struc- 
tural properties of the network 1 1], such as clustering coeffi- 
cient 1 2], average network distance |3], connectivity distribu- 
tion |4], assortativity |5], and weight distribution |6, 7, 8]. The 
effects of these properties on synchronization has particularly 
attracted the attention of researchers, partly because of the el- 
egant analysis due to Pecora and Carroll which allows us 
to isolate the contribution of tiie network structure in terms of 
the eigenvalues of the coupling matrix. 

Synchronizability of complex networks of oscillators gen- 
erally has been shown to improve as the average network dis- 
tance decreases, with one notable exception: in random scale- 
free networks, which are characterized by a strong hetero- 
geneity of the connectivity distribution |4|, synchronization 
was shown to become more difficult as the heterogeneity in- 
creases 1 10], even though the average network distance de- 
creases at the same time. Motivated by this counter-intuitive 
effect, researchers have pursued ways to enhance the synchro- 
nizability of scale-free networks by introducing directionality 
and weight to each link in the network iI^ItL IiiIi . A natural 
question arising in this context is: Given a network of oscil- 
lators with a fixed topology of interactions, which assignment 
of weights and directions maximizes its synchronizability? By 
maximization, we mean that the synchronized states are stable 
for the widest possible range of the parameter representing the 
overall coupling strength. 

The study of such a question not only provides us with 
insights into the dynamics of real-world complex networks 
but also guides us in designing large artificial networks. 
Metabolic networks — the system of hundreds of intercon- 
nected biochemical reactions responsible for the biomass and 
energy production in a cell — is a prototypic example where 
the weights and directions of feasible links (metabolic fluxes) 
are adjusted to optimize fitness, which is likely to account 
for robustness of synchronized behavior against environmen- 
tal changes 1 12|. Other examples range from the enhancement 
of neuronal synchronization for a given topology of synaptic 
connections in the brain, to the design of interaction schemes 



that optimize the performance of computational tasks based 
on the synchronization of processes in computer networks 
lITsll . The adjustment of flows in power grids and commu- 
nication patterns in social organizations are additional exam- 
ples where directional and weighted patterns can be favored 
because they can better facilitate the synchronized or coordi- 
nated behavior on which the functioning of these networks is 
based. 

Here we show that the answer to the question of maximum 
synchronizability falls outside the framework of the Pecora- 
Carroll analysis, which is built on the assumption that the net- 
work dynamics can be linearly decomposed into eigenmodes, 
i.e., the coupling matrix of the network is diagonalizable. In- 
deed, we show that maximally synchronizable networks are 
always non-diagonalizable (except for the extreme configura- 
tions where a node is connected to all the others) and can be 
constructed for any given interaction topology by imposing 
that the network: (/) embeds an oriented spanning tree, (ii) 
has no directed loops, and (///) has normalized input strength 
in each node. The fact that the networks are not necessarily 
diagonalizable has been largely overlooked in the literature, 
apparently because most previous works have focused on net- 
works of symmetrically coupled oscillators, which are guar- 
anteed to be diagonalizable. However, the same does not hold 
true in general when the network is directed, as required in 
the realistic modeling of many complex systems. Here we de- 
velop a new theory that extends the Pecora-Carroll analysis to 
the case of non-diagonalizable networks. We show that in this 
case the synchronizability is still determined by the eigenval- 
ues of the coupling matrix, but the speed at which the system 
converges toward the synchronized state may be significantly 
slower This theory is a first example of going beyond the 
traditional framework for studying complex systems based on 
either decomposition into eigenmodes or some sort of super- 
position principle. 

Consider n identical oscillators whose individual dynamics 
without coupling is governed by x = F(x), x G iR'". Now 
consider the network of these oscillators coupled via an out- 
put signal function H : J?'" J?'" along a network with a 
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symmetric adjacency matrix A = (Aij) defined by Aij = 1 
if oscillators i and j i) are connected and A^j = other- 
wise. Let Wij > denote the strength of the coupling that 
oscillator i receives from j. Thus, A represents the topology 
of interactions and W — (Wij) represents the assignment of 
weights and directions. The system of equations governing 
the dynamics of the oscillator network can then be written as 
= F(xi) + o-I]"=i AijWij[ll{xj) - H(xi)] or, equiva- 
lently. 



F(x,) -a^iyH(> 



, ,n, 



(1) 



where a is the parameter controlling the overall coupling 
strength and L — {Lij ) is the coupling matrix of the directed 
weighted network, defined by Lij — —AijWij if i j and 
I^ii = ^ Tlij^i ^ij- Note that L is not necessarily symmetric 
because the network is not constrained to be undirected. 

The maximization problem considered in this paper can be 
formulated as follows. For a given topology of interactions 
between oscillators (represented by A), we want to find the 
assignment of weights and directions (represented by W) that 
maximizes the synchronizability of the network. In order to 
address this question, we need a condition for the network 
to synchronize. For any solution x = s{t) of the individ- 
ual dynamics x ~ F(x), the completely synchronous state 
Xi = sit), i = 1, . . . , n is automatically a solution of the 
entire system Q. The question then is to determine when 
this solution is stable against small perturbations. This syn- 
chronization condition can be derived by extending the linear 
stability analysis of Pecora and Carroll |9] to the case where 
L is not necessarily diagonalizable, as follows. 

The starting point of our analysis is the observation that, 
for any n x n matrix L, there exists an invertible matrix P of 
generalized eigenvectors of L which transforms L into Jordan 
canonical form as P^^LP ~ J, where 
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and A is one of the (possibly complex) eigenvalues of L. The 
stability of the synchronous solution of Eq. is determined 
by the variational equation ^ = L'F(s)^ — (7i)H(s)^L-^, 
where ^ — (^j^, . . . , and is the perturbation to the ith 



oscillator. By applying the change of variable -q 
get 

7) = DY{s)t] - (tDU{s)t]J'^. 
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Each block of the Jordan canonical form corresponds to a sub- 
set of equations in Q. For example, if block Bj is k x k, then 
it takes the form 

171 = [DF{s) - aDU{s)]ri, (4) 
r,2 - [Z?F(s) - ai?H(s)]r,2 - aZ?H(s)r,i (5) 

r,, = [DF{s) - ai5H(s)]r,, - aDH(s)r,,_i , (6) 



where a = a\ and ry^, 773, . • . , »7fc are perturbation modes in 
the generalized eigenspace of eigenvalue A. 

For ct regarded as a complex parameter, Eq. (0} is a mas- 
ter stability equation and its largest Lyapunov exponent A(a), 
called master stability function |9|, determines the stability of 
Eq. 0: it is linearly stable iff A(ctA) < 0. The condition 
for Eq. (|5} to be stable is apparently more involved but can be 
formulated as follows. The linear stability of Eq. (0} implies 
that f]^ converges to zero exponentially as i ^ 00. Assuming 
that the norm of Z?H(s) is bounded, we have that the sec- 
ond term in Eq. (|5} is exponentially small. Then, the same 
condition A((tA) < 0, now applied to Eq. (|5}, guarantees the 
stabilizing effect of both the first and second terms, resulting 
in exponential convergence of rjj to zero as t ^ 00. The same 
argument applied repeatedly shows that 773, .... 77 j, must also 
converge to zero if A((tA) < 0. This shows that A(crA) < 
is a necessary and sufficient condition for the linear stability 
of the equations corresponding to each full block Bj. This 
condition is valid not only in diagonalizable |9] but also in 
non-diagonalizable networks. 

However, it is worthwhile noting a crucial difference be- 
tween the diagonalizable and non-diagonalizable cases. If L 
is diagonalizable, then all Jordan blocks are 1 x 1, so there 
would be no equations like Q or (|6j, and each mode of per- 
turbation is decoupled from the others. Thus, the exponential 
convergence occurs independently and simultaneously. On 
the other hand, if L is not diagonalizable, some modes of 
perturbation may suffer from a long transient. For instance, 
if we have a network of linearly coupled phase oscillators, 
di = uj — a J2j ^ij^j^ £ '5'^, then we can explicitly solve 
Eqs. 0-(|6j for the solution s{t) — uit to obtain the last per- 
turbation mode rjk = e~"* J2i=o where the constants Ci 
depend on the initial condition. Therefore, the larger the size 
k of the Jordan block, the longer the transient. 

Turning our attention back to the maximization problem, 
we first note that the eigenvalues Ai, . . . . A„ of matrix L can 
be ordered such that = Ai < Re A2 < ... < Re A„, where 
one eigenvalue is always zero because L has zero row sum 
and all the others are guaranteed to have nonnegativereal parts 
because of the Gerschgorin Circle Theorem. Thus, taking all 
the Jordan blocks into account, it follows from our stability 
analysis that the synchronous solution is stable if and only if 



A(crA,) < Ofori 2, 



(7) 



A((tAi) = A(0) > is the largest Lyapunov exponent of the 
individual oscillators and corresponds to the stability along 
the synchronization manifold. We next note that Re A2 > 
if and only if the network embeds an oriented spanning tree, 
i.e., there is a node from which all other nodes can be reached 
by followingdirected links. This condition follows from the 
recent Ref. lll4ll and generalizes the notion of connecteness to 
directed networks. We assume this condition here to ensure 
that the network is compatible with synchronization. 

In most of the previously studied cases, the master sta- 
bility function A(a), determined by F, H, and s, has been 
found to be negative in a single convex bounded region of the 
complex plane llisll . This implies the existence of a single 
interval (cTmin, Cmax) of the overall coupling strength a for 
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which synchronization is stable. Thus, the synchronizability 
of the network can be measured in terms of the relative inter- 
val (Tmax/cmin: the nctwork becomes more synchronizable as 
fmax/fmin bccomcs larger In the special case of undirected 
networks, the eigenvalues of L are real, and this measure of 
synchronizability is proportional to the ratio \2/\n O- 

A critical observation is that in order for the ratio 
Cmax/ f min to achicvc absolutc maximum for any given A(q;) 
with a convex stability region, all nonzero eigenvalues must 
be real and equal to each other The condition that the eigen- 
values must be real follows from the convexity of the stability 
region and the fact that complex eigenvalues appear in conju- 
gate pairs, while the condition that they must be equal follows 
from the fact that, for real eigenvalues, the ratio Uniax/ Cmin is 
proportional to A2 / A„. Thus, a network with 

= Ai < A2 = • ■ • = A„ (8) 

has the widest possible range of coupling strength in which 
synchronization is stable, independently of the individual 
node dynamics F, output function H, and synchronous state 
s, as long as the stability region is convex 

Under the mild assumption that the interaction topology 
allows no oscillator to interact with all the other oscillators, 
any maximally synchronizable network is necessarily non- 
diagonalizable. This comes from the fact that if L is diagonal- 
izable and satisfies the optimality condition (|8} with nonzero 
eigenvalues equal to A > 0, then all the rows of the charac- 
teristic matrix L — XI must be equal. In terms of the network 
structure, this means that each node must either have uniform 
output to all the other nodes (at least one of them must do 
so) or have no output at all. These exceptional cases include 
globally connected networks and directed star configurations. 
However, it is uncommon in a large complex network that 
an oscillator can communicate with all the other oscillators. 
Therefore, our extension of the master stability analysis to 
non-diagonalizable networks was indeed necessary to prop- 
erly address the optimization problem. 

Having observed that optimal networks are rarely diagonal- 
izable, we now show that, for any connected topology of in- 
teractions, there are assignments of directions and weights for 
which the resulting network is non-diagonalizable and max- 
imally synchronizable . We first note that maximum synchro- 
nizability can always be achieved by imposing that the net- 
work (/) embeds an oriented spanning tree, (//) has no directed 
loops, and (///) has normalized input strengths in each node, 
i.e., the total input is the same for all nodes that receive input. 
Condition (/) guarantees that Re A2 > 0, condition (//) guar- 
antees that the eigenvalues are real, and condition (Hi) then 
implies the identity (|8ji among the nonzero eigenvalues. In 
such optimal networks, we can always rank the nodes so that 
each node receives inputs only from nodes that are higher in 
the ranking (see Fig.^a) for an example). In this hierarchi- 
cal structure, information flows only from top to bottom of 
the ranking, without feedback. The optimality can be for- 
mally confirmed by noting that indexing nodes according to 
the ranking makes L a lower triangular matrix with 0, A, . . . , A 
on the diagonal, which means that A2 = • • • = A„ = A, where 
A > is the total input strength in n — 1 of the nodes. An 




FIG. 1 : (Color online) (a) Example of optimal assignment of weights 
and directions within a given interaction topology. The total in- 
put strength in each node is normalized to A, where thick, medium, 
and thin arrows indicate weight A, 2 A/3, and A/3, respectively, and 
dashed lines have zero weight. Nodes are numbered and colored to 
show the hierarchical structure in which connections are only from a 
higher level to a lower level, with no feedback loops, (b) Example of 
oriented spanning tree within the same interaction topology as in (a), 
constructed by the breadth-first search. 



important class of such maximally synchronizable networks 
consists of the oriented spanning trees themselves, where the 
normalization condition leads to uniform weights for all links 
of the tree (see Fig.^b) for an example). This example shows 
that any interaction topology admits at least n — 1, but usu- 
ally many more, optimal non-diagonalizable networks. In- 
deed, from the Matrix-Tree Theorem it follows that the num- 
ber of all oriented spanning trees is n"^2Mi' where /i2 , . • ■ , /in 
are the nonzero Laplacian eigenvalues of the underlying undi- 
rected network defined by matrix A. For a globally connected 
network, for example, the number is rt"^^, which is huge 
even for relatively small networks. All these oriented span- 
ning trees are non-diagonalizable, except for the star configu- 
ration. Oriented spanning trees can be explicitly constructed 
by the well-known procedure called the breadth-first search, 
which spans all nodes starting from an arbitrary root node. 

Physically, the optimality conditions (i)-(iii) can be under- 
stood as follows. The top node in the ranking receives no 
input and acts as a master oscillator that dominates the net- 
work dynamics. If the coupling strength a is chosen so that 
A((tA) < 0, then the oscillators that are immediately lower 
in the hierarchy and receive input from the master will syn- 
chronize themselves with the master Any oscillator receiving 
input only from these oscillators and the master must also syn- 
chronize, since normalization of the total input strength makes 
the equation effectively look as if it were receiving input from 
a single oscillator that is synchronized with the master Re- 
peating the same argument for the rest of the network, we see 
that under conditions (/)-(///) all oscillators must eventually 
synchronize and they do so for the entire range of a where 
A(crA) < 0. 

Interestingly, undirected tree networks have been found to 
be among the most difficult to synchronize |18|, in striking 
contrast to our result that directed spanning trees lead to the 
most synchronizable configurations. This highlights the sig- 
nificance of directionality of the interactions in determining 
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the synchronizability of the networks fl^ . On the other hand, 
the choice of the master oscillator in a maximally synchroniz- 
able network is completely arbitrary, despite the intuition that 
the nodes with largest connectivity would be the most natural 
choice. Moreover, the directions of the links in such a net- 
work are not necessarily related to the properties of the nodes 
they connect, even though there has been a suggestion that it 
would be related to the age of the nodes |20]. In contrast, un- 
der the stricter constraint that all feasible input connections 
have the same strength in each node, it was found f3, ^ that 
maximum synchronizability is achieved when the individual 
input strength is inversely proportional to the connectivity of 
the node, which is consistent with our result that normaliza- 
tion is key to ensuring optimality. 

The optimality conditions suggest that in designing 

a network for which synchronization is desired, it is generally 
advantageous to avoid feedback loops and to normalize input 
strength. Because these conditions typically lead to assigning 
nonzero weights only to a subset of all possible links, this in- 
teresting result can be interpreted as a synchronization version 
of the paradox of Braess for traffic flow [21], in which remov- 
ing links leads counter-intuitively to improved performance 
of the network. Furthermore, such assignment of weights not 
only maximize the synchronizability, but also minimize the 
coupling cost. The coupling cost can be defined as the sum of 



the input strengths of all nodes at the synchronization thresh- 
old |6|. If A(a) < in {ai,a2), then the coupling cost 
for any network can only be as small as ai{n — 1), which 
can be achieved by networks with global uniform coupling. 
A surprising fact, however, is that this minimum can also be 
achieved by the maximally synchronizable networks as well. 
In other words, our optimality conditions allow a network con- 
strained by an arbitrary topology to synchronize with the best 
possible efficiency. It is interesting to point out that efficiency 
optimization of traffic flow on a transportation network model 
leads to a hierarchical structure similar to that possessed by 
our maximally synchronizable networks |22|. 

Our characterization of the maximally synchronizable net- 
works can be used to test the widely assumed hypothesis that 
synchronizability plays an important role in the evolution of 
many real-world complex networks. The loop structure of 
the metabolic network of E. coli suggests that having fewer 
loops may have been beneficial for the cell (the details will be 
published elsewhere), while recent experimental findings 1 23|] 
suggest the significance of hierarchical structures in neuronal 
networks. Exploring more real data to systematically test this 
hypothesis is of critical fundamental for a better understand- 
ing of complex networks. 
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